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I. INTRODUCTION 


The time series water wave model is generally referred to 
as the Boussinesq equation. There are many versions of the 
Boussinesq equation, both the continuity equation or the 
water surface equation and the momentum equation. 
Among researchers who developed the Boussinesq 
equationareBoussinesq (1871), Dingermans(1997), Hamm, 
Madsen, and Peregrin(1993), Johnson (1997), Kirby 
(2003), Peregrine (1967), Peregrine (1972), and many 
more. These equations are formulated using continuity 
equation and Euler’sthe momentum equation, which are 
well known, but in different explanationsmaking each 
researcher has his own Boussinesq equation. 


The Boussinesq equations in this study were formulated 
using a different method from the previous Boussinesq 
equation formulation. This study formulated the continuity 
equation and momentum equation using the weighted 
Taylor series (Hutahaean, 2021a). The total acceleration in 
the momentum equation is formulated using the same 
control volume as in the continuity equation, where the 
horizontal velocity of x-direction only changes on the x- 
direction and the vertical velocity of z-direction changes 
only on the Z-axis of vertical. 
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Abstract— In this study, a Boussinesq equation was developed to model 
water waves. The equations were formulated using a weighted Taylor 
series, the Taylor series with weighted coefficients in the terms. The 
equation consists of the water surface elevation equation, the water 
particle velocity in the direction of the horizontal axis and vertical axis. 
The formulation of the continuity equation and the momentum equation is 
carried out on the same control volume with the same characteristics of 
particle velocity change. Thus, both equations have velocity variables with 
the same characteristics. The study obtained the Boussinesq equation that 
can be used with large wave amplitudes that are in line with those found in 
nature and can model shoaling and breaking well. 


The basic form of the continuity equation obtained is the 
sum of the total acceleration in the horizontal-xdirection 
with the total acceleration in the vertical Z-axis. 
Furthermore, by working on very small time intervals, the 
term acceleration versus time (localacceleration) is 
omitted, obtaining a continuity equation in the form of a 
spatial differential only, as is the case with the continuity 
equation that has been widely used. However, there is a 
weighting coefficient on the horizontal differential term. 


Furthermore, by reducing the spatial size of the control 
volume in the horizontal x-direction and in the vertical- 
zdirection, the relation between the local acceleration of 
the particle in the vertical z-axis and the local acceleration 
of the particle in the horizontal x-direction is obtained. 
This relation is used to calculate the particle velocity in the 
vertical z-direction. While the particle velocity equation in 
the horizontal x-directionis obtained from the momentum 
equation in the horizontal x-direction. The hydrodynamic 
pressure and dynamic driving force in the horizontal axis- 
xare obtained by integrating the vertical z-direction 
momentum equation on the vertical axis. 


By using the same control volume in the formulation of the 
continuity equation with the formulation of the momentum 
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equation, the definitive equation equal to the velocity of 
the particlewas obtained. Hence, there is a harmony 
between the continuity equation and the momentum 
equation. 


Continuity Equation Formulation 


The continuity equation is formulated by working on the 
principle of conservation of mass in an incompressible 
flow using the control volume in Figure (1). The input- 
output due to fluid flow in the control volume is as 
follows. 

Input, 

I = pudz+ pw 6x 

Output, 

O = p(u + du)dz + p(w + dw)dx 

As a result of the input and output, at the small time 


interval dt, there is a change in the fluid mass at the 
control volume of, 


ôm = (I — O)ôt 
j w + ôw 
u u +ôu 
ôz —> —> 


ôx 


Fig. 1: Input-output on the controlvolume 


Substitute the input and output equations, and divide both 
sides of the equation byôt ôx 6z, 
ôm ôu ow 
Stox dz Sx Gz 
For a constant control volume, the mass change in the 
control volume is 
ôm = dp 6x 6z 
equation obtained: 


6p _ ôu ôw 
st Px P oz 





; : 6 
for incompressible flow = 0, 


ôu ôw 


ox moe sees (1) 
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Equation (1) is a continuity equation in a very basic form. 
The derivative form of (1) is determined by the definitions 
of du and ow. To get the form of du and dw, the weighted 
Taylor series (Hutahaean, 2021a) is used for a function: 
f(x,z,t) = coshk(h + z) cos kx cos ot 

k is the wave number, h is the water depth, ø is the 
angular frequency, x is the horizontal axis, Z is the vertical 
axis, and t is the time. This equation is a form of velocity 
potential solution to Laplace’s equation (Dean, 
1991).cosot is used instead of sinot. The weighted 
Taylor series from the equation is: 


f(x + 6x,z + 6z,t + ôt) = f(x, z,t) 


(2) 
yis the weighting coefficient. In (2), there is a relation, 
OZ =VOxX oases. (3) 
For the horizontal velocityu(x,z,t), the only thing that 
change is the horizontal x-direction, 
u(x + 6x,z,t + ôt) = u(x, z,t) 
ty? bt + y ox E 
y ot ae ax 
The first term of right-hand side is moved to the left and 
defined, 
du = y?8t + y 6x=....44 
u=y at Y xz ) 
For the vertical velocity w(x, z, t), which changes in the 


vertical z-direction only, 


w(x,z + 6z,t + ôt) = w(x,z,t) 


ty? bt + 62 T 
EOP’ ae" ag 
dw is 
=y2 6 Ww aw 
ôw = y^ ôt ae POF ap) 


Substitute (4) and (5) to (1), 





2 5p u au" pete dw 
y“ ôt art Y ôx ax 7 ôt q t OF as 
ox 6z 


Substitutedzby (3)in the denominator of the second term 
and multiplying the equation byy 6x, 


gcc + y? bx ae + 72 ot + bz = 0 
ram Y a TSS ae 82 ae = 


dz 
...(6) 


Then, with constantéxand6z, Stis reduced to close to zero, 
the equation becomes, 
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du dw 
y? öx + ôz 750 
Z 


Substitutedzby (3)and divide the equation by y 6x, 

yet Dat) 

(7) is theweighted continuity equationthat will be used to 
formulate the water surface elevation equation, in the 
following section. 


I. THE FORMULATION OF WATER SURFACE 
EQUATIONY (x, t). 

The continuity equation (7) is multiplied by dzand 

integrated over water depth, 


1 du n 
rÍ edz + | dw =0 
-n Ax -h 


The integration of the 1* term is accomplished by the 

Leibniz integration rule. The equation of the Leibniz 

integral rule (Protter, a Morrey, &Charles, 1985)is: 
Paf , _ pP 

J PA =z u dz — fE 


Then the integration of the 1* term is, 


f #4 _ a P an dh 
Y ndx a= Vax Pa á Yin Fy YU-n Fy 


a(-h) ah 
where =-—,u 
dx dx n 


surface, whileu_pis the horizontal velocity on the sea bed. 





is the horizontal velocity at the 


This integration is substituted into(8), whilew,particle 
velocity of z-direction is substituted by weighted 
Kinematic Free Surface Boundary Condition (Hutahaean, 
2021a), which is: 


Wy yo tuys 
So, 

d rt dn, dh 
Va YUn ay Yun y 
w u 2 Wp = 0 

dt “a 


Integration is completed by working on the concept of 
depth average velocity, while w_,is substituted with the 
bottom kinematic boundary condition, 





dn df UH m dh 

Ya +’ ax +1 = yaya ax += y)aunU 5 
=0 

..(9) 


Uis the horizontal depth average velocity, „is the 
coefficient integration, His the total water depth H = h + 
n, Whilea,,is the transformation coefficient from the 
horizontal surface velocityu,to the depth average velocity 


U, whileæ„pis the transformation coefficient fromu_,to 
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depth average velocity U.These coefficients will be 
discussed in section 4. 


IW. EQUATION OF PARTICLE VELOCITY IN THE 


HORIZONTAL DIRECTION AND VERTICAL 
DIRECTION 
3.1. Momentum equation 


For the conformity with the continuity equation, the 
momentum equation will be formulated using the control 
volume in figure (1) where the horizontal velocity u only 
changes on the horizontal x-axis, while the vertical 
velocity w only changes on the vertical z-axis. Equation 
(4) is divided by dt and very small dt and dxclose to zero 
were taken, 


Du 2 du du 


=y at YU; ao (10) 





This ees is the total acceleration in the horizontal x- 
direction. Equation (5) is divided by ôt and very 
small dtandézclose to zero were taken, 

ye ey coscad (11) 


Dt at az 


This equation is the total acceleration in the vertical z- 
direction. 


A method is used to formulate the Euler’s momentum 
equationin (10) and (11),that the mass multiplied by the 
acceleration is the same as the driving force, the equations 
for the horizontal x-directionand vertical z-direction are 
obtained, respectively, 


2 du du 1dp 





Y a tY Ua pa SE (12) 
2 dw — _idp _ 
Y twos can J (13) 


ae Hydrodynamic Pressure Equation 
To obtain the pressure equationp, (13) is multiplied by 
dzand cae about the vertical Z-axis, 


rls ~ ets z| aww =-2] dp — of dz 


The integration of the 2™term on the left and the 1* and 2"! 
terms on the right is solved.The surface dynamic boundary 
conditions are donep, = Oand the equation is written as 
the equation for p, 
p dw 
== y a ar dz+= 5 (wy wy — ww) + 9(9 - 2) 
Zz 

(14) 


This equation is the hydrodynamic pressure equation. 
Next, (14) is differentiated about the horizontal x-axis, 





1 dp 2d T dw 1 dwyWy 1 daww dn 


= + + 
pax Y ax), at “T2 ax 2ax Ia 
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Hutahaean (2021b) solves the integration and differential 
on the Ist term on the right side of the equation by using 
the velocity potential of the solution to the Laplace 
equation obtaining 





dp. (= a 1 dw,w, n 1dww dn 


= = + 
pax at at) 2 dx 2 dx ax 


This equation is the equation for the hydrodynamic driving 
force in the horizontal x-direction. 


3.3. Velocity Equation in the Horizontal x-direction 
(15) is substituted to (12), 


,du yäuu „(du du 
Patra ) 
1ldw,w, 1dww dn 
2 dx 2 dx ax 
The 1‘term on the left side and the the 2™term on the right 
side cancel each other out, then this equation is multiplied 








by dz and integrated over water depth, 








du y (" duu 1dw,w. 
20%n en 
r ra. ax “TT? dx 
n dww 


+2 dz- gH ou... (16) 


-h ax 


The integration of the 2™term on the left and right sides is 
solved by the Leibniz Integral rule, 


1 duu da ø” dn dh 
| ae” = ae) dZ — Unun = Fx — Uhuh ax 


The integration of the right-hand side is solved by using 
the concept of depth average velocity, 
1 duu dfu UUH dn dh 
= dz = — Uy Uy, 7 
-n €x dx 
Puuis the coefficient integration, that will be discussed in 
section 4, then the horizontal velocity at the 
surfaceu, andthe horizontal velocity on the sea bed u_,are 


transformed intohorizontal depth average velocity, 





1 duu duu UUH dn 
K Ax dz = Ax — Ayn Aun UU mz 
dh 
Gn CunUU m> 


Qy,anda,,are transformation coefficients thatwill be 


discussed in section 4. 


In the same way, it will be obtained, 
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1 dww d WWH d 
| dz= Pow = Aww WW a 
-n dx dx dx 
dh 
Ka dx 


Pwwcoefficient integration and ay,is the transformation 


coefficient. Kinematic sea bed boundary condition 
iSW_p = —U_p z So 
dh dhy? 
WRN a RE un (a) 


ah. ah. 
als the small bottom slope, sow_p,W_p axis Very small and 


can be ignored, so 


uņin (16) is transformed into Uandw, is transformed into 
Wand divided by H, 











dU y fl duu yn @wn EWW 
Yan ae tan] a Z5- a 
1 ¢"” dww dn 
+ 2H Í ax Jax 


i i : Í a 
The result of integration is substituted tie s = dz 


and J”, = dz 
2y aU y au UUH 
Y unt 2H ax 





= Yun Cun dn _ VYCunGun UU dh _ 
2H dx 2H dx 
Lyn Gym AWW 1 AByyWWH 





2 ae ay dx 
-mm wy 2 g.(17) 


2H 


(17) is the equation for theparticle velocity in the 
horizontal x-direction. 


3.4. Velocity Equation vertical z-direction 

In (6),with constantdt,the size of6xanddzare reduced to a 
point, so that the 2"4 and 4 terms become zero, then 
divided byy? ôt, then the equation becomes, 

du rn dw 
Yat” at 

This equation is written as the equation for multiplied 
bydzand integrated about the vertical z-direction.The 


integration is completed by Leibniz integral rule and the 
A of depth average velocity is obtaining 


a 
wH + (Bw - wn )W -i = 
eTA e e oases (18) 
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In this equation is obtained from (17) while—is obtained 
from (9). The equation(18) can be said to be an expression 


of the law of conservation of energy, where there is an 
au aw an 


interaction of —, —and—. 
at’ dt at 
IV. COEFFICIENT INTEGRATION AND 
TRANSFORMATION COEFFICIENT 


In the previous sections, the integration of water depth is 
done using the concept of depth average velocity, where 
the depth average velocity (Dean, 1991) is, 


1 p 
U(x,t) =F al nude SEEN (19) 
Pu =Coefficient of integration of horizontal velocity u, 
andH =h+yn 


Furthermore, in this study, it is defined that the depth 
average velocity is the velocity at a depth of z = Zp (Fig. 
2), and the integration coefficient p4 is, 


udz 


1 n 
Bu = Hu(x, Zo t) L 









U = u(x, Zo t) 


Sea bed 


Fig 2. Depth average velocity 


Using the particle velocity equation from the velocity 
potential, 


u(x, z,t) = Gk cosh k(h + z) cos kx sin ot 
So, 
u  coshk(h+z) 


U cosh k(h + Zo) 
1 


Bu = H cosh k(h + Zo) 


n 
Í cosh k(h + z) dz 
-h 


Integration is completed and by takingy = so 


sinh k(n+4) , 
Pu = kH cosh kntz 0) 
Law of conservation of wave number (Hutahaean. 2021), 
A 
ak(h+Ż 
Aae 0...20) 
dx 
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wherek (h + 5) = kH = constant = 0m, so 


sinh 07 


i= —.....(22) 


~ @ncoshk(h+Z9)” 
So to get „is to set the value ofz beforehand. 


Ois referred to as the deep water coefficient which is 

greater than or equal to 1 where tanh@z ~1. This 

coefficient for defining the deep water limit, for 

exampleSPM (1984), using@ = 1, in this case, deep water 
0 


depth ishy = =. In this study,@ = 1.8 is used, this value 


was obtained based on a study of the solution generated by 
the model. 


n 
Í w dz = ByWH 
-h 


1 


n 
Bw = Takara, k(h + Z) dz 


By completing the integration obtaining, 


B _  Coshĝôrn-1 
W @nsinhk(h+zo)" 


..(23) 

The other integration coefficients are, 
n 

Í UU dz = B,,UUH 
-h 


1 


n 
= ———— h? k(h d 
Puu H cosh? k(h + Zo) f os i 


By completing the integration obtaining, 


_ sinh 267+67) zi 
Puu = 2 O@mcosh2 MOTRA, ) 
n 
| WW dz = B,,,WWH 
-h 


1 
~ Hsinh?k(h+z 





n 
Bww | sinh? k(h + z) dz 
0) Jen 


By completing the integration obtaining, 


> sinh 20n-On 


Bww = 


~ 26m sinh? k(h+z) ` 


(25) 


The integration coefficients are calculated in deep water. 
The value of the integration coefficient is constant, 
considering (21). 
In the momentum equation and the continuity equation, 
there are surface velocities u,,andw,, and seabed velocities 
u_p,andw__p. These variables must be transformed into 
depth average velocity of U andW, with the following 
transformation coefficients. 

cosh 61 
un = cosh k(h + Zo) 


Un = Ayn U 
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cosh 6m 


Cun = ————....(26) 


cosh k(h+zo) 
U_p = Ayn U 


1 


Qun = cosh ita 2) 
Wy = ty W 

sinh 07 
Own = sinh rmz (23) 


The vertical velocity transformation on the sea bed cannot 
be used for potential flow equations because it will 
produce a vertical bottom velocity which is equal to zero. 
For this reason, the vertical bottomvelocity is changed to 
horizontal bottom velocity by using the kinematic bottom 
boundary condition, which is, 


dh 


se =EN Ty 


dh. 
P bottom slope, so 
x 


dh 
W_p = —Ay,U ae 29) 


V. MODEL RESULTS 


First, the model is carried out on a channel with a constant 
depth, with a water depth ofh = 15.0m and a canal 
length of 300 m. The input is a sinusoidal wave with 
period T = 8 second and wave amplitude A = 1.20 m. 
The calculation constants used are y = 1.6, 0 = 1.8, and 
Zo = —0.4h. The model is executed up to 11 times the 
wave period. The model outputs are presented in Figure 3 
and Figure 4. 


In the output of the model, the resulting solution is stable 
at the execution of 11 times the wave period. This shows 
the stability of the equation and the numerical method. The 
wavelength of the model is quite short, around 25 m 
(Figure 3 and Figure 4). The particle velocity in the 
vertical direction W is greater than the particle velocity in 
the horizontal direction U with the opposite phase. The 
horizontal phase velocity U is the same as the water 
surface movement phase 7 while the particle phase 
velocity in the vertical W-direction is opposite to both. 
This shows that the change in kinetic energy in the vertical 
direction is proportional to the amount of change in 
horizontal kinetic energy with changes in potential energy 
at the water surfaceelevation. 


In addition, when U = 0 and 7 = 0, W is also equal to 
zero. This makes it easier to define the initial condition 
where att = 0, U=ņ =W 





= 0.Calculation of the water 
particlevelocity in the vertical direction using the 
Kinematic Free surface Boundary Condition, there will be 
difficulties in defining the initial value, considering that 
this equation produces a phase where at U = Oand7 = 0, 
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the particle velocity in the vertical W-direction is in the 
maximum phase. 


ANT 


0 100 200 300 


x (m) 





eta U Ww 





Fig 3. Output model on theflat bottom. 


1.5 
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£ 0.5 
= 0 
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= -1 
E 
~ -1.5 
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0 20 40 60 80 100 
x (m) 
eta U WwW 








Fig 4. Output model on theflat bottom. 


Next, the model is executed on a channel with a bottom 
slope of -0.07. The upstream water depth is 15.0 m, while 
the downstream water depth is 1.0 m, with a channel 
length of 200 m. The incoming wave is a sinusoidal wave 
with a wave period of T = 8 sec, with an initial wave 
amplitude of Ay = 1.2 m. The model output is presented in 
Fig 5. 


p 


AARAA IN 


n (m), U & W (m/sec.) 


-1 
“9 
30 50 100 150 200 
x (m) 
eta U WwW 








Fig 5. Output model on the sloping bottom. 
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Fig 5 shows that the model can model the phenomenon of 
shoaling well. At a distance of x = 150.0 m, at a water 
depth of h = 4.50 m, the wave begins to break, followed 
by the main breaking at h = 3.80 m with a breaking 
amplitude of A, = 1.5. Furthermore, after breaking, the 
model stops after one wave period. From these results, it 
can be said that the model can simulate shoaling with 
breaking even though it is not complete until the waves 
disappear. 


VI. CONCLUSION 


The formulation of the total acceleration equation in the 
momentum equation using the same control volume as in 
the continuity equation provides certainty that the velocity 
characteristics in the momentum equation and the 
continuity equation are the same. Thus, it can be said that 
the velocity stated in the continuity equation is exactly the 
generated velocity by the momentum equation creating a 
good couple between the continuity equation and the 
momentum equation. 


By defining that the depth average velocity is a 
representative velocity of the velocity at a certain depth, 
the integration coefficient, and the transformation 
coefficient can be calculated using the velocity potential 
theory. 


The weighted coefficient on the Taylor series functions to 
adjust the particlewavelength and velocity, be it horizontal 
or vertical. The greater the value of the weighted 
coefficient, the shorter the wavelength, the smaller the 
particle velocity in the horizontal direction while the 
particle velocity in the vertical direction will be greater. 


Determination of the deepwater depth represented by the 
deepwater coefficient, in addition to determining the 
stability of the solution, also determines the depth at which 
breaking begins to occur. 


In general, the model gives good results where the model 
can be carried out at large wave amplitude that match 
those in nature, and the model can simulate shoaling and 
breaking well. However, it requires development so that 
the model can simulate breaking completely with the wave 
height that gradually gets smaller until it disappears after 
breaking. 
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